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Burial of the polar magnetic field of an accreting neutron 
star. II. Hydromagnetic stability of axisymmetric equilibria 
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ABSTRACT 

The theory of polar magnetic burial in accreting neutron stars predicts that a moun- 
tain of accreted material accumulates at the magnetic poles of the star, and that, as 
the mountain spreads equatorward, it is confined by, and compresses, the equatorial 
magnetic field. Here, we extend previous, axisymmetric, Grad-Shafranov calculations 
of the hydromagnetic structure of a magnetic mountain up to accreted masses as high 
as Ma = 6 X 10~'*AfQ, by importing the output from previous calculations (which were 
limited by numerical problems and the formation of closed bubbles to A/a < IO^^A/q) 
into the time-dependent, ideal-magnetohydrodynamic code ZEUS-3D and loading ad- 
ditional mass onto the star dynamically. The rise of buoyant magnetic bubbles through 
the accreted layer is observed in these experiments. We also investigate the stability of 
the resulting hydromagnetic equilibria by perturbing them in ZEUS-3D. Surprisingly, 
it is observed that the equilibria are marg inally stable for all < 6 x IQ-^Mq] 
the mountain oscillates persistently when perturbed, in a combination of Alfven and 
acoustic modes, without appreciable damping or growth, and is therefore not disrupted 
(apart from a transient Parker instability initially, which expels < 1% of the mass and 
magnetic flux). 
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1 INTRODUCTION 

The magnetic dipole moments ^ of neutron stars are ob- 
served to decrease with accreted mass. A/a. Evidence of 
this trend is found in a variety of systems, e.g. low- and 
high-mass X-ray binaries, and binary radio pulsars with 
white-dwarf and neutron-star companions (Taam & van den 
Heuvel, 1986; Shibazaki et al., 1989; van den Heuvel & 
Bitzaraki, 1995) , although there is some debate over whether 
the trend is monotonia (Wijers, 1997). Several theoreti- 
cal mechanisms have been proposed to explain the /i(Ma) 
data, including accelerated Ohmic decay (Urpin & Gep- 
pert, 1995; Urpin & Konenkov, 1997), fluxoid- vortex in- 
teractions (Muslimov & Tsygan, 1985; Srinivasan et al., 
1990), and magnetic screening or burial (Bisnovatyi-Kogan 
& Romberg, 1974; Romani, 1990). With regard to the lat- 
ter mechanism, Payne & Melatos (2004) (hereafter PM04) 
calculated a sequence of two-dimensional, hydromagnetic 
(Grad-Shafranov) equilibria describing the structure of the 
magnetically confined mountain of material accreted at the 
magnetic poles of the neutron star. The mountain is con- 
fined by magnetic stresses near the equator, where the field 
is compressed (Melatos & Phinney, 2001). These solutions 
are the first of their kind to explicitly disallow cross-field 



transport of material as the mountain evolves from its ini- 
tial to its final state (cf. Mouschovias, 1974), as required in 
the ideal-magnetohydrodynamic (ideal-MHD) limit. PM04 
found that is screened substantially above a critical ac- 
creted mass A/c ~ 1O~^M0, well above previous estimates 
of Afc < 10"^" A/0 (Hameury et al., 1983; Brown & Bildsten, 
1998; Litwin et al., 2001). 

PM04 calculated equilibria up to A/a < IO'^A/q, falling 
short of the mass required (~ O.IA/q) to spin up a neu- 
tron star to millisecond periods (Burderi et al., 1999). This 
is supplied by a large class of mass donors like LMXBs 
(Strohmayer & Bildsten, 2006), even given nonconservative 
mass transfer (Tauris et al., 2000). Grad-Shafranov calcu- 
lations are stymied above A/a ~ 1O~'*M0 by physical ef- 
fects (e.g. magnetic bubbles form above the stellar surface) 
and numerical effects (e.g. steep magnetic gradients hin- 
der iterative convergence). In this paper, we extend the 
^(Ma) relationship to Ma ~ 1O~^M0 by loading equilibria 
with A/a ~ lO"*A/0 into ZEUS-3D, a multi -purpose time- 
dependent, ideal-MHD code for astrophysical fiuid dynam- 
ics, and adding extra mass quasistatically through the outer 
boundary of the simulation volume. 

PM04 also left open the important question of the sta- 
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bility of the hydromagnetic equilibria. Distorted magnetic 
configurations, in which the polar flux is buried beneath the 
accreted overburden and compressed into a narrow belt at 
the magnetic equator, are expected prima facie to be un- 
stable. Indeed, the Grad-Shafranov analysis in PM04 hints 
at the existence of an instability by predicting (informally) 
the existence of magnetic bubbles as steady-state solutions. 
In this paper, we systematically explore the stability of the 
equilibria by evolving them in ZEUS-3D, subject to linear 
and nonlinear perturbations. 

The structure of the paper is as follows. In section 2, 
the necessary theory is summarised and the solution method 
described. The formalism of PM04 is again used here. The 
numerical Grad-Shafranov solver is described in appendix B 
of PM04 and appendix C of Mouschovias (1974). In section 
3, we verify ZEUS-3D against a set of test cases relevant 
to the problem of magnetic burial; the implementation is 
described thoroughly in Appendix A. In section 4, we explore 
the late stages of magnetic burial (10^^ < M^/Mq < 10"^) 
by adding mass quasistatically to equilibria from PM04 in 
ZEUS-3D. In section 5, we discuss the linear and nonlinear 
stability of the equilibria in the regime 10~® < Ma/M© < 



10" 



The paper concludes, in section 6, with a discussion 



of the limitations of our analysis and suggestions for future 
numerical work. 



and Rt is its radius, mass-flux conservation in ideal MHD 
provides the integral constraint 

X |27r j dsrsin6l|VV'r^exp[-(0- </>o)/c^] 



(2) 

We prescribe the mass-flux distribution to be dM/dtj) = 
(Ma/2V'a) exp(— i/'/V'a), where Va = ip*R*/Ra. is the flux 
enclosed by the inner edge of the accretion disk at a distance 
-Ra and -i/i* is the hemispheric flux. For the boundary condi- 
tions, we fix xp to be dipolar at r = i?, , assume north-south 
symmetry, fix the tp = field line, and leave the field free at 
large r. 

Equations (1) and (2) are solved numerically using an 
iterative relaxation scheme (Mouschovias, 1974; Payne & 
Melatos, 2004). The mean residual as a function of iteration 
number is shown in figure 1(b), corresponding to the Grad- 
Shafranov equilibrium for m = Ma/Mc = 0.16 displayed in 
figure 1(a). The form of F{tp), found from (2), varies from 
Ma = to Ma = 0.16Mc in the manner displayed in figure 
1(c). 



2 PHYSICS OF MAGNETIC BURIAL 
2.1 Axisymmetric equilibria 

During accretion onto a neutron star from a binary com- 
panion, matter piles up on the polar cap, funnelled by the 
magnetic tension of the polar magnetic flux tube. Once Ma 
exceeds ~ 10~^Mo, the hydrostatic pressure at the base of 
the accretion column overcomes the magnetic tension and 
matter spreads over the stellar surface towards the equator, 
dragging along frozen- in polar field lines (PM04). The dis- 
torted magnetic field leads to screening currents which act 
to decrease the magnetic dipole moment. Figure 1(a) illus- 
trates the magnetic 'tutu' formed for Ma = W~^Mq, cut 
off at ten density scale heights. The polar mountain of ac- 
creted material (dashed contours) and the pinched, fiaring, 
equatorial magnetic belt (solid contours) are plainly seen. 

At first glance, one might expect such equilibria, with 
their steep density and magnetic field gradients, to be un- 
stable. Interestingly, this expectation is largely false, as we 
show in section 5. We summarise the key equations and no- 
tation of magnetic burial here (PM04). 

The steady-state ideal-MHD equations for an isother- 
mal atmosphere (p = c^p, i.e. adiabatic index 7 = 1) reduce 
to the force equation Vp + pV0 — /io~^(V x B) x B = 0. 
For an axisymmetric configuration in spherical polar coordi- 
nates (r, 6,(j)), a flux function tp('r, 9) generates the magnetic 
field B via B = V?/'(r, 6')/(r sin 61) x e^. The flux function 
satisfies the Grad-Shafranov equation 

A"-0 = F'(V') exp[-(<^ - 4,0) li], (1) 

where p, p, 4>, 4>o, Cs, and denote the matter density, 
pressure, gravitational potential, surface gravitational po- 
tential, sound speed, and Grad-Shafranov operator respec- 
tively (PM04). In the limit ho = cIrI/{GM») < R», 
<p = GMtv/R^, where M« is the mass of the neutron star 



2.2 Critical accreted mass Mc 

According to naive estimates based on hydromagnetic force 
balance between ^q"^(V x B) x B and Vp at the polar cap, 
the magnetic field (and hence p,) is distorted appreciably 
away from its initial configuration for Ma > Mc ~ 10~^Mq 
(Brown & Bildsten, 1998; Litwin et al., 2001). However, self- 
consistent solutions of (1) and (2), in which mass does not 
migrate across fiux surfaces and the back reaction from equa- 
torial magnetic stresses is included (PM04), predict a larger 
value of Mc, given by 

^ = ■nGM.BlRl/{2p.octMQ) 

Mq 

- — '^(lifc) {-Ik) 

(ro) (io6i^s-i) • 

In the regime Ma < Mc, the Green function analysis in 
PM04 gives 

V(r,e) = Vd(r,e)[l - mb\l - e"^/"")] , (4) 

with m — Ma/Mc and 6 = 4'*/'<Pa.- The distorted field lines 
develop a large tangential component as mass is added, such 
that |B| increases substantially for A/a ^ M^/a, with a = 
R*/ho^ The self-consistent density distribution associated 
with (4) is given by 

p{r,e) = poF[i,{r,e)]e-''^''° (5) 



^ The extra factor of a comes from differentiating ?/' with respect 

to e. 
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Y'/Y', Scale heights 

Figure 1. Hydromagnetic equilibrium for a = 50 (scaled from real stellar dimensions as discussed in section 3.3) and Ma/Mc = 0.16, 
with Gx = Gy = 64, showing (a) magnetic field lines (solid) and density contours (dashed) and (b) the corresponding mean residual 
versus iteration number for an underrelaxation parameter = 0.99. (c) Final F(ip) (solid) and initial F(^p) (dashed) versus V- (d) 
Magnetic dipole moment (normalized by its surface value) as a function of altitude above the stellar surface (in units of ho)- 



with 

^^^^ = 2^ exp(-^)(l - i,/b)'/'[J(i;)r\ (6) 

where wo write F — Fhf)Cs /M^, pa = A4^/hf), and ip = ip/ips^, 
and the form factor satisfies J(V') ~ 1 to better than 10 
per cent for all 6 except near the equator, 89.5° < < 
90° (see figure Al of PM04). In the small-Ma hmit, F is 
calculated by assuming ij} to be dipolar and evaluating the 
(small) correction from (1) using Green functions. 



where /ii = ijjfR* is the initial dipole moment. Equation (8) 
depends only on m = M^jM^, not Ma and Mc individu- 
ally, just as in (4). Figure 1(d) shows how ii decreases as 
a function of altitude due to the screening currents within 
the first few density scale heights. It is very important to 
note that deviates substantially from its dipole form for 
Ma/Mc > 6~^, whereas /i deviates from /ii for Ma/Mc > 1, 
which occurs at a much later stage of accretion (because 
6^1 usually) . 



2.3 Screening the magnetic dipole 

The magnetic dipole moment is defined in terms of the radial 
component of the magnetic field by 

H=^r^ J d(cose) coseBr{r,e), (7) 

assuming axisymmetry. In the regime Ma < Mc from (4), 
we obtain 

H/Hi = (1 - Ma/Mc) , (8) 



3 NUMERICAL SIMULATIONS OF 
MAGNETIC BURIAL 

Distorted MHD equilibria like the one pictured in figure 1(a) 
are notoriously unstable. In this section, we describe how the 
astrophysical MHD code ZEUS-3D can be used to investi- 
gate the stability of our equilibria. We discuss the set up 
of ZEUS-3D in section 3.1 and appendix A, some verifica- 
tion experiments in section 3.2 and appendix A, and the 
curvature rescaling required to render the magnetic burial 
problem tractable in section 3.3. 
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3.1 ZEUS-3D 

ZEUS-3D is a multipurpose, time-dependent, ideal-MHD 
code for astrophysical fluid dynamics which uses staggered- 
mesh finite difi^erencing and operator spHtting in three di- 
mensions (Stone & Norman, 1992a). In this paper, we re- 
strict the dynamics to two dimensions, disabling the third, 
but employ a spherical polar grid, appropriate for an ax- 
isymmetric stellar atmosphere. Previous numerical work fo- 
cused on the magnetic poles of the accreting star (Brown 
& Bildsten, 1998; Gumming ct al., 2001); here, by con- 
trast, equatorial magnetic stresses arc treated fully by sim- 
ulating a complete hemisphere. The density and magnetic 
field strengths are read into ZEUS-3D from the output of 
our Grad-Shafranov code (PM04). The time-step Atz is 
set by the Gourant condition satisfied by the fastest MHD 
wave modes. Details regarding the parameters, initial con- 
ditions, boundary conditions, verification tests, and coordi- 
nate choices in the runs are given in appendix A. 



3.2 Verification 

Before implementing the burial problem in ZEUS-3D, we 
ran a sequence of simpler verification cases. 

First, we reproduced the classical results for the nonlin- 
ear evolution of the Parker instability of a plane-parallel field 
in rectangular geometry (Mouschovias, 1974). We achieved 
agreement on the minimum stable wavelength Acrit ~ 
4:nho[B^ / (4np) + 1]^^'^'^ and reproduced the final, nonlin- 
ear ly evolved equilibrium state to an accuracy of 5 per cent. 

Second, to test spherical coordinates in ZEUS-3D, we 
studied the evolution of a spherical isothermal atmosphere 
containing zero magnetic field, followed by an atmosphere 
containing a dipolar magnetic field; both these states are 
analytic force-free equilibria [(V x B) x B = 0)]. ZEUS-3D 
confirmed that these are indeed equilibria; they do not al- 
ter significantly even after several thousand Alfven times. 
The condition at the outer boundary r = rm is chosen to 
suit the problem at hand. The outflow condition leaves the 
magnetic fleld free but allows some mass loss, which we min- 
imize by keeping rm large, to prevent the atmosphere from 
evaporating over long time-scales. The inflow condition arti- 
ficially pins the magnetic field, thereby introducing a radial 
magnetic field at the outer boundary which increases arti- 
ficially when integrated at rm- It is used at an intermediate 
stage in the bootstrap algorithm (described in section 4.2) 
when adding mass in the Ma 2> Mc regime. 
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Figure 2. Magnetic dipolc moment, fi, as a function of accreted 
mass, m = M^/Mc, for b = Rn/Ra, = 3 and a = Ra./ho = 
1.86 X 10* (crosses), a = 50 (triangles), a = 100 (pluses), a = 
500 (diamonds), all found with the Grad-Shafrajiov code, and 
a = 50 with mass added through the outer boundary in ZEUS- 
3D (squares). 



fixed. This has the advantage that the minimum density in- 
creases as a decreases (because the atmosphere extends fur- 
ther), decreasing the Alfven speed and hence increasing the 
ZEUS-3D time step Atz- We set a — 50, a good compromise 
that ensures a 3> 1 (as for a realistic star) while keeping the 
computational burden reasonable. This choice corresponds 
to a model star with = lO'^M© and = 27 m. Our 
numerical results confirm that /j, is independent of a (see 
section 4.1), as predicted analytically. This is illustrated in 
figure 2, where /i(Ma) is plotted for several a values; even 
for Ma < Mc, the deviations are less than ten per cent. 

Our ZEUS-3D grid reaches an altitude Xm = rm — R* = 
Who (cf. > Vfho in PM04). We adopt a restricted domain 
for two reasons: (i) to maximize grid resolution near the 
surface of the star, where gradients are steepest; and (ii) 
to stop the time-step, Atz, from becoming so small that 
run time and numerical dissipation become excessive (as 
discussed in appendix A). According to the Gourant con- 
dition, Atz scales as wXinax pmin oc e^^"'''''' . Note that 
the pole-to-equator sound and Alfven crossing times (aho / Cs 
and aho/vA respectively) also decrease as a decreases. 



3.3 Curvature rescaling 

The characteristic radial (ho) and latitudinal (R*) length 
scales are very different in a neutron star, creating numerical 
difficulties which must be handled by rescaling the problem. 
For a typical neutron star with Cs = 10^ m s~^ (Brown & 
Bildsten, 1998), we find ho = 0.54 m, a = R*/ho = 1.9 x 10"*, 
and To = ho/cs = 5.4 x 10~^ s, where to is the sound crossing 
time over a hydrostatic scale height. 

Two input parameters are varied in the simulations: Ma 
and b. In the regime Ma < Mc, where (4) holds, ip and con- 
sequently n depend only on m = Ma,/Mc and not explicitly 
on Ma, suggesting that we can artificially reduce a = R»/ho 
by reducing M* and R*, as long as we keep ho oc tit/M^ 



4 LATE STAGES OF MAGNETIC BURIAL 

(Ma » Mc) 

In this section, we investigate the evolution in time of the 
highly distorted Grad-Shafranov equilibria that arise for 
Ma ^ Mc ~ 1O~*M0. We aim to answer a question on 
which the Grad-Shafranov analysis is silent: what happens 
when significantly more than W^'^Mr.) is accreted at a rate 
that is slow compared with the Alfven time and instability 
oscillation period (see section 5)? To do this, we employ the 
bootstrap approach described in section 4.2. 
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4.1 jj, versus Ma 

Our numerical results confirm that the magnetic dipole mo- 
ment is essentially independent of a — R,/ho. Figure 2 dis- 
plays jJL{Ma) for several a values, with Ma > Mc achieved 
by adding mass through the outer boundary of the simula- 
tion box (squares). As predicted analytically, there is negli- 
gible dependence on a for Ma < Mc. For Ma > Mc, devia- 
tions of less than a factor of two occur. The outer boundary 
condition (radial B) artificially increases integrated at 
Xui = lOfto, by about 10 per cent; this is considered further 
below [see figure 3(f), where /i is plotted as a function of 
altitude x (dotted curve) for m = 0.16]. Hence the magnetic 
dipole moment plotted in figure 2 is an upper bound. 



Errors are < 1% when compared against other runs 
with Sm > lO/io- For example, for m = 0.16 and Gx = Gy = 
128, the minimum dipole moment is 0.8695 at Sm = 10 and 
0.8555 at im = 20. (A 64 x 64 grid gives similar results, e.g. 
H = 0.875 at = 10.) 

4.2 Bootstrap accretion: Ma » Mc 

The data in figure 2 extend to Ma « 5Mc, yet the Grad- 
Shafranov code only produces equilibria for Afa ;S 1O~*M0, 
and breaks down when magnetic bubbles first appear at 
Ma > Mc&~^ (section 4.3). This state of affairs is un- 
satisfactory because, in many real accreting systems (e.g. 
LMXBs), Ma exceeds substantially the critical mass for 
bubbles to form (Taam & van den Heuvel, 1986; van den 
Heuvel & Bitzaraki, 1995). For such systems, the previous 
calculations are useful for the early stages of accretion, at 
Ta = Ma/Ma < 10''6~^ yr, but teach us little about the final 
stages of accretion and hence the relic magnetic structure 
once accretion stops. 

In this section, we perform a numerical experiment to 
address this issue. We begin with a hydromagnetic equi- 
librium configuration for p and where Ma is just below 
the critical value for bubbles, calculated using the previous 
Grad-Shafranov method (PM04). We load this configura- 
tion into ZEUS-3D. We then add mass on the polar flux 
tube < V V'a, at a rate slower than all the hydro- 
magnetic time-scales involved in the process of coming to 
equilibrium (including the Alfven time-scale and the period 
of the global Parker oscillations discussed in section 5), but 
at a rate faster than the true accretion rate Ma (so as to 
make the computation tractable). We then see what final 
magnetic structure we get after a realistic amount of mass 
is accreted, e.g. Ma ~ O.IM© for LMXBs. 

This approach does not track properly any processes 
that operate on time-scales between ta and Ta, such as Hall 
drift (Gumming et al., 2004) and Ohmic difi^usion (Romani, 
1990). However, it docs model all timo-dcpcndcnt ideal- 
MHD effects in the context of polar magnetic burial for 
the first time. Note that this experiment would be com- 
pletely impractical without the Grad-Shafranov equilibria 
previously computed for Ma6 < 1O~*M0 (PM04), because 
the separation of the hydromagnetic and accretion time- 
scales is too great if one starts from Ma = 0, even when 
the accretion is accelerated artificially. 

Figure 4 illustrates one bootstrapping cycle; it is also 



a good illustration of the dynamics of magnetic burial for 
Ma » Mc. Starting from a steady state equilibrium with 
m = 1.6, we add mass to the outer boundary, with inflow 
boundary conditions. The magnetic field lines are tied to the 
outer boundary by this condition in ZEUS-3D, preventing 
H from decreasing there and raising artificially closer to 
the star. To overcome this, we change the outer boundary 
condition from inflow (B pinned) to outflow (B free) after 
adding progressively more mass through the polar flux tube, 
with Ma(6') = Ma.maxe"'"'"'' ^ and Ma.max ^ 2x IO^^Mo/to 
when scaled up from the model (o = 50) to a realistic neu- 
tron star. Mass falls along field lines [figure 4(a)] towards 
the stellar surface and flattens the equatorial 'tutu', until 
the magnetic tension matches the hydrostatic pressure. Mat- 
ter subsequently fiows oquatorward, dragging magnetic field 
lines with it [figure 4(c)]. Superposed on this process are 
sound and Alfven oscillations, which are discussed in sec- 
tion 5; they cause wiggles in the field lines, visible above the 
spreading matter (2 - 3 m above the surface). The oscilla- 
tions show up clearly in the time evolution of the magnetic 
dipole moment, and mass quadrupole moment, plotted in 
figure 5. [FYirther discussion of the mass distribution ap- 
pears in Mclatos & Payne (2005).] Finally, the mass inflow 
is stopped and the outer boundary condition is switched 
to outflow, allowing the magnetic field lines to relax [figure 
4(c)]. The configuration in figure 4(e) becomes the initial 
state for further mass to be added. 



The dipole moment attained after adding a total of 3 x 
lO^^M© on a 128 x 128 grid is shown in figure 6 as a function 
of altitude. We find that a minimum dipole moment of w 
3 X 10~^ times the surface value is obtained for m = 4.76 
at an altitude ~ ^o- Smaller dipole moments arc expected 
for larger Ma, but numerical constraints prevented us from 
improving the grid resolution nearer the equator sufficiently 
to see this. 

One might be tempted to ask whether the intermediate 
step of relaxing the field at fm is really necessary. It is, and 
figures 3 and 4 illustrate why. If too much mass is added, 
it is found numerically that the field lines break off the un- 
derlying magnetic dipole and, when allowed to relax, remain 
pinned to the outer boundary. Figure 5 shows as a func- 
tion of time, measured ^-^ 1.5 m above the stellar surface 
where the spreading mostly occurs, jj, reaches a minimum 
after ~ IOOto, even when more mass is added, because of 
line-tying at rm- Adding more mass is ineffective because 
closed magnetic loops form. They are apparent in figure 7, 
which shows the configuration at f = 400 (starting from the 
m = 1.6 equilibrium). This configuration cannot be used as 
the initial state for the next bootstrapping iteration because 
the magnetic field lines remain fixed at the outer boundary. 



To validate the bootstrapping method, we take an equi- 
librium that is free of magnetic bubbles (i.e. m <^ 1) and 
add mass up to a value of Ma for which the Grad-Shafranov 
equilibrium is already known. Generally speaking, we repro- 




Figure 3. Magnetic field evolution with mass added to the outer boundary, (a), (c), (e): Magnetic field lines {solid) and density contours 
{dashed) for t = 0, 500, 2000, corresponding to m = 0.16,0.22,0.4 respectively, (b) Dipole moment as a function of time, measured at 3 
m above the stellar surface, (d) Quadrupole moment as a function of time, (f) Dipole moment as a function of height for t = {dotted), 
500 {dashed), and 2000 {solid). 



duce the Grad-Shafranov value of to an accuracy of 5 per 
cent. 

Figure 3 shows the results of such a test, beginning 
with m = 0.16 and adding mass at a rate Ma.{0) = 
Ma.maxe-*"^'"' ^ with Ma.max w 7 X IQ-^Mq/tq. We point 



out several important features, (i) The initial n in figure 
3(f) rises towards the outer boundary (see section 3.2). (ii) 

fi is unchanged at the outer boundary, because the mag- 
netic field lines arc tied there by the flux-freezing condition, 
to accomodate the inflowing mass. The magnetic fleld lines 
bend towards the equator, and n decreases closer to the stel- 
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Figure 4. Illustration of one step in the bootstrapping algorithm. Magnetic field lines (solid) and density contours (dashed) (left) and 
the magnetic dipole moment as a function of height (right) for m = 1.6 at t = (top), m = 2.2 at t = 100 (middle), reached by adding 
mass at the outer boundary, and m = 4.1 at t = 400 (bottom). 



lar surface, (iii) Closer to the surface, /i reaches a minimum 
within 5 per cent of that obtained from the Grad-Shafranov 
equiUbrium. (iv) The kinks apparent in the magnetic field 
lines are a result of numerical dissipation on the grid scale. 



4.3 Bubbles 

When the Grad-Shafranov equation is solved with ip free 
(Neumann) at the outer boundary, closed bubbles of mag- 
netic field, disconnected from the inner (Dirichlet) bound- 
ary, can arise (PM04) for Ma > l.BMcfe-^ Physically, they 
are generated when the toroidal screening current exceeds 
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Figure 5. Evolution of the magnetic dipole moment, as measured 1.5 m above the stellar surface {left), and the mass ellipticity (right), 
when mass is added to a m = 0.16 equilibrium. 




Figure 7. (Left.) Magnetic field lines (solid) and density contours (dashed) when mass is added to a m = 0.16 equilibrium, reaching 
m = 4.1 at t = 400. (Right.) Magnetic dipole moment as a function of altitude when the m = 4.1 configuration at t = (dotted) relaxes 
in the presence of an outflow outer boundary condition, at t = 500 (dashed) and t = 900 (solid). 



a threshold. However, the field hnes are deformed continu- 
ously from a simply connected initial state (e.g. a dipole, or 
something close to it is the best guess for a typical neutron 
star) in ideal MHD. Therefore, the bubbles are not realizable 
in an accreting neutron star even though they are admissible 
mathematical solutions to the steady-state boundary-value 
problem. Instead, the bubbles point to a loss of equilibrium, 
analogous to that which occurs during eruptive solar phe- 
nomena (Klimchuk & Sturrock, 1989), where no simply con- 
nected hydromagnetic equiUbrium exists. In the language of 
the Grad-Shafranov analysis in section 2, the source term 
oc F'(ip) in (1) increases with Afa&, boosting A'^ip and cre- 
ating flux surfaces with ?/; < or i/; > which cannot 
connect to the star and either form closed loops or are an- 
chored at infinity (here, the accretion disk). 

How docs a bubble evolve in ZEUS-3D? We import a 
Grad-Shafranov equiUbrium for Ma = 1.6Mc into ZEUS-3D, 
retaining the self-consistent density p = F[ij;{r,6)]e~'^^'^° in 
the region < ip < ip^ but replacing it with an isothermal 



atmosphere inside the bubble (ip < and tp > tpt), where 
strict flux-freezing would imply p = (matter cannot enter 
without crossing flux surfaces). Figure 8 shows part of the 
evolution for m = 1.6. The bubble rises to rm, during the 
first Alfven oscillation of the magnetosphere, at the Alfven 
speed. During the bubble's rise, 0.5 per cent of the accreted 
mass (4.4 x 10^^ Mq when converted to reahstic neutron star 
parameters; see appendix A6) is ejected through the outer 
boundary in one Alfven time. This compares with thermal 
evaporation of less than 0.001 per cent in the same time, 
which occurs naturaUy given outflow boundary conditions 
(see appendix A). For m > 1.6, the bubble oscillates about 
an equilibrium point rather than rising buoyantly. 



4.4 Uniform density increase 

Another route to computing equilibria with Ma > lO^^M© 
is to start from a bubble-free Grad-Shafranov equilibrium 
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with Ma < Mc and increase the density uniformly across 
the grid by some factor (usually ^ 1), while leaving the 
magnetic field unchanged. If we start from the initial dipole 
(Ma = 0), this kind of numerical experiment is badly 
controlled; it leads to excessive numerical dissipation and 
mass loss through r = Vm- However, if we start from 
A/a ~ 10~^Mq, the readjustment is gentler. During the 
early stages of the readjustment, some classic instances of 
the Parker instability are observed. For example, figure 9 
shows how the initial state for m = 1.12 evolves after the 
density is increased five-fold uniformly to m = 5.6. The mag- 
netic field component Be is compressed in the r direction on 
an Alfven time-scale — ripe conditions for the Parker insta- 
bility. The blistering becomes clear after 160 Alfven times 
(third frame). After 280 Alfven times, the material settles 
down to a new equilibrium. Importantly, negligible mass and 
magnetic flux (less than 1 % of the total) is lost through this 
blistering. 

The magnetic dipole moment //. plotted in figure 10, 
settles down to « 0.06/ii. This is larger than we expect for 
m — 5.6, given that we found /i(m = 4.7) < 0.01 using 
the bootstrapping method. Two factors explain this: (i) the 
sudden density increase does not allow the magnetic field 
and matter to gradually find their equilibria, and (ii) flux- 



freezing is violated by increasing the density as above. This 
method can be used in conjunction with bootstrapping with 
less severe density increases (factor <C 5). 



4.5 Size of the polar cap 

In our numerical experiments, we choose the polar cap ra- 
dius to be relatively large because (i) it helps the Grad- 
Shafranov solution converge numerically while still captur- 
ing the essential idea of a polar mountain, (ii) it raises the 
minimum density pmin oc e"** on the grid, preventing Atz 
from becoming too small via the Courant condition, and 
(iii) in a real star, the polar plasma flow can spread due to 
Rayleigh- Taylor and Kelvin-Hclmholtz instabilities (Arons 
et al., 1984). We have not managed to get ZEUS-3D run- 
ning for b ^ 3. For 6 = 3 and 6 = 5, /i reaches an equilibrium 
value in a few hundred Alfven times, and the oscillation fre- 
quencies are roughly the same, with — > 0.85, 0.82 for 
6 = 3, 5 as t — » inf. 
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Figure 9. Evolution of tlie Parker instability for Ma = 1.12Mc with the density increased uniformly across the grid by a factor of 5. 
Snapshots of the magnetic field lines are shown at t = 0, 120, 160, 200, 240, and 280 Alfven times (top left to bottom right). 



4.6 Equatorial magnetic belt 

Melatos & Phinnoy (2001) predicted that the magnetic field 
at the equator intensifies during magnetic burial, if mag- 
netic flux is conserved. The results presented here confirm 
this. As seen in figure 1(a), the magnetic field is "combed" 
away from the pole, and flattened against the stellar sur- 
face towards the equator, increasing Bg at the expense of 
Br- The maximum magnetic fleld strength, Bmax, is given 
as a function of Ma in Payne & Melatos (2006b) Note that 
the analytic approximation in the limit of small Ma captures 



the equatorial belt provided that J{i)) is calculated properly 
(PM04). 



5 STABILITY 

In this section, we investigate whether the Grad-Shafranov 
equilibria of magnetically conflned mountains calculated in 
this paper (Ma > Mc) and PM04 (Ma < Mc) are stable 
to small and large perturbations. To do this, we load the 
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Figure 6. Dipolc moment versus altitude when mass is added 
to the configuration in figure 4(c), aX t = (dotted), t = 100 
(dashed), and t = 500 (solid), corresponding to m = 2.23,2.86, 
and 4.76 respectively. 
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Figure 10. Magnetic dipolc moment as a function of time during 
the blistering of the Parker instability displayed in figure 9. 

numerical output from our Grad-Shafranov code into ZEUS- 

3D, evolve it forward in time, and report on the nature of 
any instabilities observed. In ideal MHD, any instabilities 
manifest themselves as slow, Alfvcn, or fast magnetosonic 
waves, modified by buoyancy effects in a stratified medium. 

We explore the stability of the system in the context 
of three numerical experiments: (i) we test how the small 
perturbation evolves, which arises from the numerical er- 
ror in converting from the Grad-Shafranov to the ZEUS-3D 
grid (section 5.1, 5.2 and 5.4); (ii) we consider the fate of 
large perturbations (section 5.5); and (iii) we compare the 
evolution in ZEUS-3D with the convergence of the Grad- 
Shafranov algorithm. 

5.1 Linear stability: numerical experiment 

Equilibria of the kind depicted in figure 1(a), where the mag- 
netic field is markedly distorted, are normally expected to 



be unstable. To assess this here, we track separately the evo- 
lution of Wg, VFk, and Wb, the gravitational, kinetic, and 
magnetic energies respectively, as the magnetic mountain 
evolves in ZEUS-3D. The total energy is given by 

W = WB + Wg + Wk, 

where we define 

Wg = J d^^p4>, Wk = J d^x(pvV2), 

and 

Wb = J d'x(BV2Mo). 

Key observables of the system, such as the magnetic dipole 
moment fi, defined in (7), and the mass cUipticity e, axe 
also tracked. We do not plot the total thermal energy 
Wp = f d^xPlogP in what follows because ZEUS-3D holds 
it constant during isothermal (but not adiabatic) runs. 

We begin, as an example, with the equilibrium for 
rn = 1.6 [figure 4(a)], our starting point for adding mass 
through the outer boundary in section 4.2. Figure 11 dis- 
plays Wg, Wk and Wb a function of t for this case. The 
energies are normalized by their maximum values so that 
the energy exchanges are clear, because their absolute val- 
ues differ by several orders of magnitude: Wg dominates, 
with Wb = 3.7 x W^Wg and Wk = 4.2 x lO'^Wg ini- 
tially. After 300ta, we obtain Wb = 3.4 x 10"^ Wg and 
Wk = 1.2 X 10"'"^Wg. 

The evolution proceeds as follows. The equilibrium state 
imported into ZEUS-3D is not exact due to grid resolution, 
imperfect convergence in the Grad-Shafranov code, and nu- 
merical discrepancies when translating between the Grad- 
Shafranov and ZEUS-3D grids (see appendix A). Initially, 
during the first oscillation, Wg and Wb are converted to 
Wk. When the oscillation overshoots, the energy fiow re- 
verses direction. In ideal MHD, where there is no dissipa- 
tion, the oscillation persists. In ZEUS-3D, where there is 
some numerical dissipation and energy is radiated by hydro- 
magnetic waves through the outer boundary, the oscillations 
are damped, eis in figure 11. However, this damping is slow. 
(Experiments demonstrating this numerical dissipation are 
reported in appendix A.) Note that Wk oscillates at twice 
the frequency of Wg and Wb, because the leading term in 
Wk is quadratic in perturbed quantities (oc 6v^) whereas Wg 
(cx 5p) and Wb (oc B5B) are linear. 

The dipole moment, normalized to its value at the stel- 
lar surface, is plotted in figure 12 for m = 1.6. The oscil- 
lation period equals that of Wb for the reason above. The 
kinks in the curve occur when the magnetic field reflects off 
the boundaries at the pole and equator. The oscillations are 
identified as acoustic modes (see section 5.2). 



5.2 Global MHD oscillations 

We explore here the physical nature of the oscillations re- 
sulting from small perturbations to the hydromagnetic equi- 
libria and how their amplitude and period depend on Ma. 

Figure 13 shows the time evolution of n and mass el- 
lipticity for several Ma values and logarithmic scaling of the 
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Figure 11. Kinetic (solid), magnetic (dashed), and gravitational 
potential (dotted) energies, normalized to their maximum values 
as a function of time for m = 1.6, Gx = Gy = 128. 
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Figure 12. Magnetic dipole moment normalized to its surface 
value, plotted as a function of time for m = 1.6, Gx = Gy = 128. 

altitude. For reference and numerical comparison, the re- 
sults for a linear scaling in altitude are shown in figure 14. 

The key result is that the equilibria are marginally stable: 
the buried field is not disrupted significantly, but the con- 
figuration oscillates about its equilibrium state. Two modes 
are clearly present: (i) a short period oscillation, with fixed 
period 13ro for all Ma (and thus p), which is an acoustic 
mode with velocity Cg; and (ii) a longer period oscillation, 
with period increasing with Ma as displayed in figure 15, 
which is an Alfven mode. The Alfven oscillation frequency 
is fitted by 

/a « 0.003To-'(Ma/Mc)-''"Hz (9) 

for a = 50, which, when scaled to a realistic neutron star, 
yields 

/a « 17(Ma/Mc)-'/'Hz. (10) 

There are a few numerical considerations. The magnetic 
dipole moment artificially rises from its equilibrium value 
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Figure 13. .Vlagiictic dipole moment (top) and mass quadrupolc 
moment (bottom) as a function of time for M^/Mc = 0.05.3, 0.16, 
0.32, 0.48, 0.64, 0.8, 0.96, 1.12, 1.6, 2.4, 3.2, and 4.0 (top to 
bottom for /j, and bottom to top for the quadrupolc moment). 
Grid resolution; Gx,y = 128, with logarithmic scaling in altitude. 

towards the surface value for 0.2 < w- < 1.6. under the fol- 
lowing conditions, (i) There is insufficient resolution close to 
the stellar surface, e.g. a linear grid in r with Gx = 128 is 
not good enough. A comparison of figures 13 (logarithmic r 
grid) and 14 (linear r grid) brings out this point, /i and the 
quadrupolc moment exhibit oscillations, whose amplitude 
grows and whose mean varies for a linear grid, whereas the 
amplitude is damped and the mean remains constant for a 
logarithmic grid, (ii) The outer boundary condition is set to 
outflow, leading to a small amount of mass loss. While the 
mass loss amounts to < 2% hy t = 10^, it allows the field 
lines at the magnetic equator to flare up towards the pole 
and destabilize the equilibrium. These problems are over- 
come by (i) using a logarithmic scale in r to increase the 
grid resolution near the stellar surface, (ii) managing care- 
fully the translation from logarithmic scaling in the Grad- 
Shafranov grid to the ZEUS-3D grid (see appendix A), and 
(iii) setting the outer boundary condition to inflow. 
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Figure 14. Magnetic dipolc moment (top) and mass quadrupole 
moment {bottom) as a function of time for Ma/Afc = 0.053, 0.16, 
0.32, 0.48, 0.64, 0.96, 1.12 (top to bottom for n, and bottom 
to top for the quadrupole moment). Grid resolution: Gx,y = 64, 
with linear scaling in altitude. 



5.3 Identifying MHD modes 

To classify the oscillation modes discovered in the previous 
section, we need the dispersion relation of small-amplitude 
MHD waves in the limit (verified by the simulations) where 
the wavelength of the p and B perturbations is small com- 
pared with [p|/[Vp| and [7/>|/|V'!/'|. The force equation gov- 
erns stability: intuitively, if a fluid (Lagrangian) displace- 
ment ^ creates a force F in the opposite direction, equi- 
librium tends to be restored and the system is stable. 
Upon analysing the perturbation in Fourier modes, ^(r, t) oc 
^(k, a))e~'^'''''~"'*\ the force equation reads 



Pow^4 = -c,Vokk-^(k,a;)+ /Xq '{[k x (kx (^ x Bo))] x Bq} . 

(11) 

In the WKB limit k ^ |V-!/'|/|V'|, the dispersion relation 
(11) supports two modes: shear Alfven waves, with 



Figure 15. Oscillation period (in units of the Alfven time) as 
a function of Ma, with linear (crosses) and logarithmic (dia- 
monds) grid sampling in altitude. The Alfven period is well fit by 
300To(Ma/Mc)^/^. The period of the sound mode (triangles) is 
also plotted. Note that to is defined at a particular value of Ma. 
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The shear Alfven wave is incompressible (k • ^ = 0) . It prop- 
agates analogously to a transverse wave along a string under 
tension, the magnetic field lines playing the role of the string. 

Let us denote the right-hand side of (11) by F(^). The 
operator p~^F is linear. It can be shown [see p 244 of Goed- 
bloed & Poedts (2004) for a proof] that in ideal MHD, the 
operator p^^F is also self-adjoint, so that the eigenvalues o;^ 
in (11) are real (making uj either real or purely imaginary). 
Two classes of solution occur: (i) stable pure waves (w^ > 0), 
and (ii) exponentially growing instabilities (a;^ < 0). There 
are no damped oscillations in ideal MHD, even in a nonuni- 
form background; any dissipation observed is numerical. 

A full theoretical calculation of the two-dimensional sta- 
bility is beyond the scope of this paper. However, we at- 
tempt to interpret the numerical results displayed in equa- 
tions (9) and (10) in terms of known results for a gravitat- 
ing, magnetized slab (Goedbloed & Poedts, 2004). We as- 
sume that the slab is infinite and homogeneous in the y- and 
^-directions, contained between two planes at a; = a;i and 
X = X2; that the equilibrium p and tp vary in the a;-diroction; 
and that the excited modes satisfy the boundary conditions 
(,x{xi) = £,x{x2) ~ 0. The results for a slab provide insight 
into the oscillation modes of a neutron star's magnetosphere, 
where the above boundary conditions are well satisfied due 
to a combination of line tying at r = R, and stratification 
at all r (the component of ^ perpendicular to equipotential 
surfaces is small). Working with coordinates in a field- line 
projection, the Alfven and slow continuum frequencies are 
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estimated to be (Goedbloed & Poedts, 2004) 

2 i -1 -1 17,2 2.-1 IP -1ei2 

Wa = 47r//o P , t^s = 47r//o — rr^P F 

7P + 47r//o S2 

(15) 

with F = — iB • V. Prom our numerical simulations, we find 
empirically that « 1/2 wavelength of the dominant oscilla- 
tion mode fits within one quadrant, yielding radial and lat- 
itudinal wavenumbcrs ~ 2/i?, . This gives F ~ 2Bma,x/R* 
and thus u)% = 4B^ax/(^*Pmax), using the majcimum mag- 
netic field strength and density as characteristic values. 
Empirically, from our simulations, wo find Bmax ~ 10^ T 
for Ma = Mc, pn,ax « 10" (M,/Mc) kg m-^ and hence 
/a = u>A/{2n) = 32(Ma/Mc)"^/^ Hz. This scales the same 
way as the numerical result /a = 17{Ma./Mc)~^^^ Hz but is 
larger, due to the crudeness of our slab approximation and 
numerical estimates. The plasma varies with density (and 
Ma) from being high-/3 at the stellar surface to low-/? far 
from the star, given that Cs is uniform where the mass is 
concentrated and the waves are launched. 

Figure 15 shows the oscillation period as a function of 
m. When analyzing the MHD evolution, wo extract the pe- 
riod of the waves by Fourier analysis in order to distinguish 
which modes are excited. For the Alfven mode, as the av- 
erage density increases (with m), va oc MaT^^^^ decreases, 
and thus the oscillation period increases. The Alfven period 
(in units of the Alfven time) is fitted by 300T()(Afa/Mc)"^/^. 
The period of the acoustic mode remains constant through- 
out all the simulations we perform. 



5.4 Transient Pcirker instability 

In the experiments described in section 4.4, it is observed 
that increasing the density uniformly by a factor > 4 causes 
Parker instabilities to occur. Figure 9 consists of a sequence 
of frames illustrating the evolution of the Parker instability 
for an equilibrium state whose density is uniformly increased 
five-fold. The wavelength A of the instability subtends « 1 
radian, i.e. A « 0.3i?*. The radial density profile is exponen- 
tial, while the magnetic field is mostly confined to a layer of 
thickness ~ ho, directed perpendicular to VP. This is the 
classic situation in which Parker instabilities are expected 
(Mouschovias, 1974). 

The evolution of /i, as measured at a; = Who, is dis- 
played in figure 10. Note that, after the Parker blister sub- 
sides, fi returns to its original value. The equilibrium is not 
disrupted permanently; indeed, less than ~ 1% of the ac- 
creted layer and frozen-in magnetic flux are expelled from 
the simulation domain. In this respect, the instability can 
be considered transient. 

The Grad-Shafranov equilibria imported from PM04 
are generated by an algorithm that can follow, by succes- 
sive relaxation, the full nonlinear evolution of the Parker 
instability (Mouschovias, 1974; Parker, 1966). Therefore it 
is unsurprising that the equilibria arc stable; the relaxation 
algorithm evolves the magnetic field quasistatically through 
a sequence of intermediate "states" quite close to those that 
the real solution to the time-dependent MHD equations 
would pass through. True, the final state is distorted, and 
one might ask why the buoyancy of the compressed mag- 
netic flux does not drive long-wavelength, slow MHD modes 
that overturn the accreted matter on the Alfven time-scale 



— something that does not occur in the ZEUS-3D runs (ex- 
cept for the transient in figure 9). The explanation is that 
the equatorial magnetic 'tutu' represents the end state of 
a Parker instability that occurs quasistatically as we add 
material. The tutu is the analogue of the stable magnetic 
'blisters' which form when the plane-parallel Parker insta- 
bility saturates, while the polar mountain is the analogue of 
the material that drains into the magnetic valleys. 

To test the assertion (PM04) that the iterative solution 
of the Grad-Shafranov equation coupled to a fiux-frcczing 
condition provides a good proxy for the true time-dependent 
behaviour, we begin with a magnetic dipole with excess mass 
loaded into the polar fiux tube in both the Grad-Shafranov 
code and in ZEUS-3D. In ZEUS-3D, this is effectively the 
same as applying a uniform density increase (section 4.4) to 
the m = equilibrium, ending up with m = 1.6 in this exam- 
ple. We plot /u as a function of time (ZEUS-3D) and iteration 
number (Grad-Shafranov) in figure 17. An underrelaxation 
parameter Q = 0.99 is selected by trial and error to give a 
comparable rate of convergence in the two codes. The initial 
progress to equilibrium is similar in both cases. However, 
after 100 iterations (or equivalently t > 100), /x fluctuates 
more in ZEUS-3D than in the converging Grad-Shafranov 
code. In ZEUS-3D, the added mass initially squashes the 
fleld into a layer ~ 1/4 the thickness of the equilibrium 
layer [figure 4(a)]. The field then bounces back and transient 
Parker-like instabilities occur, as we expect (p increases by 
a factor > 4). This is best illustrated in figure 18, which 
shows a snapshot of the configurations at t = 100 (in ZEUS- 
3D) and iteration 100 (in the Grad-Shafranov code). Notice 
how the field lines are bent by the transient instabilities in 
ZEUS-3D but remain smooth in the Grad-Shafranov code. 
This comparison also illustrates why Grad-Shafranov equi- 
libria are a necessary starting point for ZEUS-3D when aim- 
ing to reach stable equilibria for large Ma > 10~*Mq. Note 
that the convergence-matching shown here is not a proof of 
the stability of the equilibria (Asseo et al., 1978). 



5.5 Large perturbations 

Next, we investigate what happens when the Grad- 
Shafranov equilibria are perturbed far from equilibrium, i.e. 
5B/B ~ 1; the perturbations considered thus far are small, 
arising mainly from imperfect translation between the grids 
of the two codes. 

We perturb the magnetic field, while simultaneously re- 
specting the boundary conditions, by setting 

B B{1 + dsm[iv(r - R*)/{r^ - R*)] sinO} . (16) 

Figure 16(a) shows the fate of this perturbation for different 
amplitudes 5. Note that the perturbation does not strictly 
respect flux-freezing because it changes dM/dip slightly. 
Nonetheless, even for significant perturbations, the equilib- 
ria are not disrupted — a significant and robust result. For 
S < 0.2, /i settles back to within 20 per cent of its initial 
value. For 5 > 0.2, /u does not settle back to its initial value, 
but it does settle down to some steady value. Again, line- 
tying plays a key role in conferring such exceptional stability 
on the system. The evolution also depends on the sign of 5. 
The perturbation is a low-order spatial mode (in r and 9) 
and either compresses {6 > 0) or relaxes {6 < 0) the mag- 
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Figure 16. {Top.) Magnetic dipolc moment as a 
function of time for perturbation amplitudes <S = 
-1.0, -0.5, -0.25, 0, 0.25, 0.5, 1.0 {bottom to top), when 
the magnetic field is perturbed as dcscrilxxi in (16) but p 
is not altered. {Bottom.) As above, but with the magnetic 
field perturbed via (17) and then iterated once through 
the Grad-Shafranov code to obtain the self-consistent p, for 
S = -0.2, -0.1, 0, 0.1, 0.25, 0.5 {bottom to top). 

netic field, with compression causing a back reaction which 
increases /u [as can be seen in figure 16]. 

The above perturbation (16) induces an efi^ective redis- 
tribution of mass in flux tubes. To avoid this, we perturb the 
magnetic potential tp in the Grad-Shafranov code according 
to 

^^^{1+6 sin[7r(r - i^,)/(rm - -R*)] smO} (17) 

and iterate once to calculate p self-consistently without 
changing dM/dtjj, so that the above mass redistribution is 
performed self-consistently. We then load the self- consistent 
perturbed equilibrium into ZEUS-3D as before. The results 
are shown in figure 16(b). For a given 5, the amplitude of the 
oscillations is reduced, confirming that the mass redistribu- 
tion contributes to the back reaction. The main effect is to 
bring the final value of ^ closer to its initial value for S > 0.5. 
The evolution stiU depends on the sign of 5, suggesting that 
mass-fiux redistribution remains somewhat important. 
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Figure 17. Magnetic dipole moment as a function of time in 
ZEUS-3D {solid) and as a function of iteration number in the 
Grad-Shafranov code with underrelaxation parameter = 0.99 
{dashed). 



6 CONCLUSIONS 

We find that, when the MHD equilibria computed in PM04 
are perturbed, they are marginally stable. The magnetic 
field remains essentially undisrupted, with the exception 
of transients that expel a tiny fraction of flux and mat- 
ter. MHD oscillations, including both acoustic and Alfven 
modes, cause /j, and e to oscillate about their mean values. 
Using ZEUS-3D, we extend the results of PM04 to larger Ma 
in two ways, (i) We add mass through the outer boundary at 
an artificially accelerated rate which is still slow compared 
to the MHD equilibration time-scale and oscillation period 
(as in a real neutron star) . We use a bootstrapping method 
in which the mass is progressively added then allowed to 
settle to equilibrium, (ii) We increase the density instanta- 
neously uniformly, and then allow the coufiguratiou to settles 
to equilibrium. The equilibria we obtain are stable to Parker 
modes, as one might expect given that they are the output 
of a Grad-Shafranov code which follows the full nonlinear 
evolution of Parker instabilities. 

The existence of stable magnetic mountains with per- 
sistent (albeit oscillating) mass quadrupole moments (e < 
10~®) raises the prospect that accreting millisecond pulsars 
(Wijnands & van dor Klis, 1998) are sources of gravitational 
radiation. This application is pursued in Melatos & Payne 
(2005) and Payne & Melatos (2006a). Likewise, the existence 
of a stable equatorial magnetic belt of intense magnetic field, 
with its ability to impede thermal transport, has interesting 
implications for the persistence of millisecond oscillations in 
the tails of type I thermonuclear X-ray bursts in LMXBs 
(Muno et al., 2002). This application is pursued in Payne & 
Melatos (2006b). 
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6.1 Limitations 

The work presented above and in PM04 can be generalised 

by relaxing certain key assumptions. Some areas to explore 
include resistive efltects (e.g. Ohmic diffusion) (e.g. Konar & 
Bhattacharya, 1997), sinking of accreted material through a 
"soft" stellar surface (e.g. Konar & Choudhuri, 2004), Hall 
currents (e.g. Rheinhardt et al., 2004), nonaxisymmetry, and 
other equations of state. For example, by assuming a hard 
surface, wc end up with densities at the base of the accreted 
layer that can be unphysically high. And it is well known 
that MHD equilibria can be stable in two dimensions, yet 
unstable in three dimensions, even with line tying. The mag- 
netization of the accreted plasma also warrants inclusion; it 
has been ignored by all authors except Uchida & Low (1981). 

The main bottleneck hindering progress towards Grad- 
Shafranov equilibria for Ma > 1O~^M0 is grid resolution at 
the magnetic equator. At present, we believe that the only 
way to alleviate this is by appropriate logarithmic gridding 
in e. 

Ohmic dissipation due to electron-phonon and electron- 
impurity scattering (e.g. Bhattacharya & Srinivasan, 1999) 
is not modelled in this paper, to keep the problem manage- 
able. Ultimately, it should be. It is important for magnetic 
structures ~ 1 m in size, which develop for Ma > 1O~^M0 
(Melatos & Payne, 2005; Brown & Bildsten, 1998; Gumming 
et al., 2001). Note, however, that nonideal effects are already 
present at some level in our calculations accidentally: the 
grid introduces numerical dissipation, because a field line is 
defined by linear interpolation between grid points and the 
error in this interpolation effectively allows matter to move 
across field lines. 



and many simulations of toroidal and nonaxisymmetric field 

structures appear in the plasma physics literature. These 
can be studied in ZEUS-3D and will be explored in a future 
paper. 

Line-tying boundary conditions, as distinct from peri- 
odic boundary conditions, change the character of the basic 
MHD waves (Gocdbloed & Halberstadt, 1994). Line-tying 
boundary conditions generally conflict with the phase rela- 
tionships between the components of the displacement vec- 
tor ^ of the three pure modes. MHD waves of a mixed 
nature occur as a result. Stability is enhanced, because it 
takes extra energy to bend the longitudinal component of 
the magnetic field. We do not compare line-tying and peri- 
odic boundary conditions quantitatively here but point out 
that line tying enhances stability even in three dimensions 
(Gocdbloed & Halberstadt, 1994). 

Short-wavelength ballooning instabilities are stabilized 
by line-tying until the overpressure at the top of the 
stellar crust exceeds the magnetic pressure by a factor 
8i?* arcsin(6"^/^)//io (Litwin et al., 2001). While they are 
predicted to occur within Iiq of the surface, we do not ob- 
serve them in these simulations. Resistive ballooning (Velli 
& Hood, 1986) and resistive Rayleigh- Taylor modes (Khan 
& Bhatia, 1993) are not considered by virtue of our re- 
striction to ideal MHD, but they are known to be com- 
mon in distorted axisymmetric equilibria similar to those 
in Figure 1(a). These instabilities can allow plasma slip- 
page faster than Ohmic diffusion, but slower than the Alfven 
time.^ To properly model these instabilities globally, large 
toroidal quantum numbers need to be considered (Huang 
et al., 2006). 



6.2 Instabilities in three dimensions 

This paper must be generalised to three dimensions before 
the stability question is truly settled. As a trivial example, 
by restricting ourselves to two-dimensional equilibria with 
= 0, we suppress unstable modes involving the toroidal 
magnetic field. The Parker instability in the galactic disc 
has been studied in three dimensions by Kim et al. (1998), 
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APPENDIX A: MAGNETIC BURIAL IN 
ZEUS-3D 

In this appendix, wo explain briefly how to simulate the 
problem of magnetic burial in ZEUS-3D, so that the reader 
can reproduce and generalize our results.^ For details on the 
general design and operation of ZEUS-3D, please consult 
Stone & Norman (1992a,b). 

After discussing the control variables in ZEUS-3D, we 
provide several test cases for reference. The Parker insta- 
bility in rectangular coordinates, whose nonlinear evolution 
was computed by Mouschovias (1974, 1996), tests the ba- 
sic functionality of ZEUS-3D. An unmagnctized, isothermal 
atmosphere in spherical coordinates tests the boundary con- 
ditions on p and v. A magnetized, isothermal atmosphere in 
spherical coordinates tests the boundary conditions on B 
and yields the equilibrium state for Ma = 0. 

Al Control variables 

The number of active grid zones in each coordinate, 
ggenl : nbl and ggen2 : nbl for i and j respectively, is Gx x 
Gy (typically 128 x 128). In addition, there are two ghost 
zones at each edge of the grid, for setting boundary condi- 
tions. The zones in, say, the i direction (r in spherical co- 
ordinates) can scale geometrically using ggenl :xlrat, such 
that the width of zone (i+1) is xlrat times zone i. This is 
how we implement logarithmic coordinates to increase grid 
resolution near the stellar surface, and to match input data 
from the Grad-Shafranov code in PM04 (section A7). There 
are several conditional compilation switches which control 
the geometry. To implement two dimensions, enable KSYM, 
and set the number of zones in the third ordinate, ggen3:nbl 
to one. For Gartesian coordinates, enable XYZ; for spherical 
polar coordinates, enable RTP. When printing out the grids, 
we use nxpx = nypx = 400 pixels and suppress the third 
dimension {z or 0) by setting pixcon: ipixdir = 3. 

The results in this paper assume an isothermal equation 
of state for simplicity, thus ISO is enabled. Self-gravity, GRAV, 
is disabled. Dimensionless quantities are chosen as follows: 
(So = 1 (automatic in ZEUS-3D), G = 1 (set grvcon:g=l), 
Cs = 1 (set eqos : ciso=l), and ho = \ {ho = c-^Ri /GM,). 
With these choices, the basic units of mass, magnetic field 
and time become Mo — hocl/G, Bo = [/ioc*/(G'/io)]^''^i and 
To = ho/cs- In ideal MHD, the evolution is controlled by 
the geometry as well as the ratio of magnetic to thermal 

^ This work also serves as a prototype of a simulation pipeline 
being developed by the Australian Virtual Observatory. We have 
demonstrated a preliminary version of this pipeline in which hy- 
dromagnetic equilibria generated by the time-independent, Grad- 
Shafranov code described in PM04 are output in a standard for- 
mat (VOTable) and fed into ZEUS-3D to be evolved in time. 
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pressure, a = B'^ / [noci p); the physical units of B and p 
do not enter separately. Varying G, set to one by default, 
changes the unit conversions but not the physics. ZEUS-3D 
may be run as a hydrodynamic code without magnetic fields 
by disabling MHD, reducing the run time. 

The time-step Atz is adjusted adaptively so that the 
Courant-Priedrichs-Lewy stability condition is satisfied in 
every zone, with the tolerance determined by the Courant 
number (hycon:courno = 0.5 by default). In the case of an 
isothermal atmosphere and subsonic fluid motions, it is set 
by max;(cs,VA). Since one has va oc 1/p ex in a typical 
isothermal atmosphere, the code is limited in the range of 
altitudes and densities that it can faithfully simulate. The 
grid is dumped at times separated by pixcon : dtpix and the 
code runs for a maximum time pcon:tlim. 

A2 Verifying the Parker instability in a 
rectangular geometry 

A test case which is easy to implement (and relevant to 
the problem of magnetic burial) is the Parker or magnetic 
Rayleigh- Taylor instability (Parker, 1966; Mouschovias, 
1974, 1996). The initial equilibrium state is a semi- infinite, 
exponentially stratified atmosphere in a uniform gravita- 
tional field, with a magnetic field perpendicular to the grav- 
itational force. This configuration is unstable to overturn 
(slow MHD) modes longer than a critical wavelength Acrit- 
We simulate the Parker instability in the region —X < 
x<X,0<y<Y. A uniform gravitational acceleration gy 
in the y direction is implemented in ZEUS-3D by placing a 
point mass far from the simulation box, at y = Rgs ^ Y. 
The boundary conditions are set to be periodic at x = ±X 
(niib, noib = 4) and reflecting a,t y = 0,Y (nijb = nojb 
= !)• 

In dimensionless units, the equilibrium initial conditions 

are 

p{x, y) = 0.5 exp[-y/(l + a)] (Al) 

and 

^(x,y) =2(l + Q)exp[-y/2(l + a)], (A2) 

with = (V X A)^ = -e-*'/2(i+"). The initial velocity per- 
turbation is taken as Vy — esm{ny/Y)cos(-Kx/X), whore, 
for e = 0.3, the kinetic energy equals 0.23% of the magnetic 
energy; cf. Mouschovias (1996). 

The stability condition is tested by importing a stable 
equilibrium state, with X < Acrit /2, and checking that the 
kinetic energy remains negligible compared to the magnetic 
and gravitational potential energies. The ratio of the kinetic 
energy to the total energy is given in table Al, scaling as R~g 
for Rgs ^ 10^. We choose .Rgfi = 7.6 x 10* for most runs. 
This test also confirms that the equilibrium is stable for 
{X,Y) — (7, 25), just below the critical value Acrit/2 = 7.26 
(Mouschovias, 1974), and is unstable for {X,Y) = (8,25). 

We also set (X,Y) — (12,25) and verify the output 
against the results given in figure 3 of Mouschovias (1996). 
We choose the same units, namely Cs = 6.2 km s~*, ho = 
1.3 X 10** m. To = ro/cs = 6.5 x 10*^ yr. Mo = 5.8 x lO" 
kg, po = 2.9 X 10"^^ kg m"^ Bo = 3.7 x 10~** T. (To 
have exactly the same units as Mouschovias (1996) choose 
G = 6.67 X 10"" in ZEUS-3D.) Figure Al shows the 



ReS Wb 



7.6 X 10^ 


-1.03 X 10^ 


12.9 


1.03 X 10" 


-1 


7.6 X 102 


-1.06 X 10* 


14.0 


1.75 X 10- 


-3 


7.6 X 10^ 


-1.06 X 10^ 


14.0 


1.32 X 10" 


-5 


7.6 X 10^ 


-1.06 X 10^ 


14.0 


1.22 X 10" 


-6 


7.6 X 10"' 


-1.06 X 10^ 


14.0 


1.69 X 10" 


-6 


7.6 X 10^ 


-1.06 X 10* 


14.0 


1.75 X 10- 


-6 



Table Al. Gravitational, magnetic and kinetic energies (dimen- 
sionless units) resulting from a point mass located at ReS such 
that M/tt^^ = 1. The initial state is the equilibrium given by 
(Al) and (A2), with f = 10 and {X,Y) = (7,25), perturbed as 
described in the text. 



ZEUS-3D results at the three times displayed in figure 3 
of Mouschovias (1996). The results are in good accord, dif- 
fering by less than 5 per cent. A mountain is clearly present 
where the magnetic field bends down in the negative y di- 
rection. 



A3 Spherical, isothermal atmosphere, zero 
magnetic field 

Having verified a known solution, we now move closer to 
our goal and consider spherical geometry (RTP). We consider 
first an isothermal atmosphere with no magnetic field to test 
the efi^ect of the boundary conditions. In spherical geometry, 
the curvature of the star, characterized by a = R*/ho ~ 
2 X 10*, enters the problem. To cope with this large value, 
we scale coordinates in ZEUS-3D logarithmically (section 
Al). However, the code halts prematurely when the time- 
step Atz becomes too small; for a > 10^, it does not run at 
all. The problem is most severe far from the surface, where p 
is small, va is large, and Atz is small. As an initial tost case 
in spherical coordinates, we sot up a spherical isothermal 
atmosphere with a negligible magnetic field (a = 4x 10~*) in 
2D (single quadrant). The initial condition is an exponential 
isothermal atmosphere 

p{r,9) = exp[— o(r — a)/r] w exp(— x) (AS) 

threaded by a dipolar magnetic vector potential 

A^(f,6>) = 10-*°f-^sin6i (A4) 

The boundary conditions are reflecting at the pole {Bg = 0, 
nijb = 1) and at the equator {Br = 0, nijb = 5), dipolar 
at the surface (niib = 3, i.e. inward flow at zero velocity, 
fixed density, dipolar magnetic field), and free at the outer 
boundary (f = fm) (niob = 2, outward fiow). Note that, in 
ZEUS-3D, a fixed magnetic field must be accompanied by 
infiow, oven if the inflow is at zero speed. 

The experiment described in this section aims to cali- 
brate the role of the outer boundary condition. Equations 
(A3) and (A4) almost describe a force-free equilibrium ex- 
cept that matter can evaporate through the outer boundary. 
However, alternative boundary conditions offered in ZEUS- 
3D are not appropriate, e.g. inflow artificially pins the mag- 
netic field lines at the boundary. The amount of mass loss 
is monitored and fm is chosen large to minimize it. There is 
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Figure Al. Nonlinear evolution of the Parker instability, showing 
magnetic flux surfaces (solid) and isodensity contours (dashed) 
for (X, Y) = (12, 25), f = 10, 20, 30, and Gx = Gy = 64 



a trade off between minimizing mass loss as well as simula- 
tion run time (determined primarily by Atz), as illustrated 
in tabic A2. Logarithmic grid scaling has little effect on the 
amount of mass which evaporates through the outer bound- 
ary. The outer boundary condition artificially increases the 
dipole moment at large r, as discussed in the main text 



and illustrated in figure 3(f). The run time is reduced if the 
magnetic calculations are switched off by disabling MHD; the 
results are indistinguishable. The results in table A2 are for 
= 6.4 X 10® m, M« = 6.0 x 10^" kg, Cb = 290 m s'^, 
ho — 8.4 X 10"^ m. With those choices, the basic units be- 
come Mo = ci/G = 1.2 X 10^^ kg, po = Mo/rl = 2.1 x 10^ 
kg m-^ Bo = [tJ.oct/{G4)]'-^^ = 15 T , TO = ro/cs = 29 
s; the alert reader will note that these parameters describe 
the Earth! 

A4 Spherical, isothermal atmosphere, dipole 
magnetic field 

As a prelude to testing the stability of the Grad-Shafranov 
equilibria found by PM04, we consider the simpler situa- 
tion of the spherical, isothermal atmosphere (section A3) 
threaded by a substantial dipolar magnetic field. Note that 
this is a valid forcc-froo equilibrium state [(V x B) x B = 
in (1) except at r = 0], except for evaporation. This config- 
uration corresponds to the pre-accretion neutron star with 
A/a = 0. It serves to calibrate ZEUS-3D and estimate the 
magnitude of numerical errors. 

The relevant neutron star parameters are Rt = 10^ m, 
M, = 1.4Mo, Cs = 10® m s~^. With these choices, the basic 
units become Mo = c^/G = 1.5 x 10^^ kg, po = Mo/rl = 
9.5 X 10^2 kg m-^ Bo = [^.oct/(Grl)f'^ = 3.5 x lO" T , 
To = ro/cs ~ 5.4 X 10^ s. However, as discussed in section 

3.3, we consider a scaled version of the real star with a — 50, 
leaving the physical scale height, ho = CsRt/(GM^) = 0.54 
m, unchanged. 

To estimate the errors in the code, we calculate the 
magnetic dipole moment and mass quadrupole moment and 
compare to the theoretical values fi = 0.5o^B* and e = 
in table A3. For b = 1, we obtain e = 1.02619; for 6 = 3, 
e = 1.23373. The errors are less than 1 per cent. Figure A2 
shows how the dipole and mass ellipticity vary with time 
for these cases. The dipole moment is accurate to < 0.2%. 
The minor variations are caused by numerical jitter (also see 
figure A5) and evaporation at the outer boundary (section 

5.4. The dependence of the error on grid size is tabulated in 
table A3. 



A5 Outer boundary 

The outflow condition at the outer boundary has minimal 
effect on the evolution if the outer boundary is sufficiently 
far away. To test this, we try several values of fm, keeping 
fm small enough so that the Alfven time-step, set by the 
Courant condition with va oc B^ /p, is not too short. 

Numerical dissipation occurs in ZEUS-3D for runs 
longer than > 10^ time-steps. Kinks in the magnetic field 
appear at grid points, as illustrated in figure A5. If the outer 
boundary is set too many scale heights away, the minimum 
density across the grid is very low, the maximum Alfven 
velocity is very high, and the time-step, Atz, determined 
by the Courant condition, is very short. This means that 
(i) the simulation takes a very long time to run, and (ii) 
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Figure A2. The dipole (left) and mass ellipticity {right) as a function of time for an isothermal atmosphere with a = 100 and rm = 110 
in a 64 X 64 grid. Parameters: po = 1 and Bo = 0-1) xlrat = 1.03. 



the numerical dissipation causes the simulation to stop. For 
a = 50 and fm = 70, the code stops at f < 2; for fm = 80, it 
stops at f < 0.2. For a = 100 and fm ~ 110, the code stops 
at f = 10.26, but with ggenl:xlrat = 1.03 instead of 1, 
it continues until f = 27.22. For o = 1000 and fm = 1010, 
the code stops immediately. Throughout the paper, we use 
fm ~ 60, so that the runs are long enough to give sufHcient 
time to assess the stability of the equilibria. Table A4 shows 
how the grid ratio affects the time to halt. 



A6 Converting to realistic a values and between 
codes 

Here we give the formulae for converting between the Grad- 
Shafranov code and ZEUS-3D. Furthermore, we detail the 
effects of converting from a small (a = 50) to a realistic 
(a w 10*) star. With ho fixed, but allowing a to vary, we 
have = aho = 27(a/50) m, M^/Mq = a?hocl/{GMQ) = 
1.0138 X 10"^(a/50)^ and 

Ms^^-'""'" IsoJ 1 10^^) • 

(A5) 



a 


fm 


% mass loss 


run time (min) 


50 


60 


0.5 


5 


50 


55 


14.1 


5 


50 


52 


17.0 


5 


20 


30 


1.72 


3 


40 


50 


1.2 


4 


100 


110 


0.3 


10 


500 


510 


0.1 


20 



Table A2. Proportion of the initial mass lost from the simulation 
volume, and run time for i = 100, as functions of the dimension- 
less hydrostatic scale height . The grid size used was 50 x 100. 
Naturally, the relative run times are more useful than the abso- 
lute. 



ho is the dimensionless unit of length. The dimension- 
less units in the Grad-Shafranov code (subscript 'G') are 
po,G = Ma//io = 9.8 X 10-^'^(m/0.16)(a/50)* kg in-'\ 
Bo,G = B*aV(2fe) = 4.2 x 10^''(a/50)^(fe/3r^ T. In ZEUS- 
3D (subscript 'Z'), po,z = ci/{Ghfi) = 9.6 x 10^^kgm-^ 
Bo,z = {no/Gh^y/^ci = 3.5 x lO" T. For conversion 
from the equilibrium code to ZEUS-3D, the factors are 
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Figure A3. The dipolc (left) and mass cUipticity (right) as a 
function of time for an isothermal atmosphere with a = 500 and 
fm = 510 in a 64 X 64 grid. Parameters; po = 1 and Bg = 0.1, 
xlrat = 1.03. 



Figure A4. The dipole (left) and mass ellipticity (I'ight) as a 
function of time for an isothermal atmosphere with a = 40 and 
fm = 50 in a 64 X 64 grid. Parameters; po = 1 and Bq = 0.1, 
xlrat = 1.03. 



grid size 




e 




64 X 64 


0.999101 


4.32268 X 10- 
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1.26503 X 10" 


-4 



Table A3. Errors as a function of grid size as measured by the 
deviation of /i and e from theoretical values for //t and et = for 
a = 100, fm = 110 and f = 0. 



PG,z = Po,g/po,z = GMa/c? = 1.3 X 10-^(m/0.16)(a/50)*, 
and Bg,z = Bo,g/Bo,z = B.a7(26c2)(G/ig/Ato)'/'' = 1-2 x 
10-*(o/50)2(6/3)-i. 

A7 Logarithmic coordinates 

PM04 concentrated maximum grid resolution near the equa- 
tor and stellar surface where Vp and Yip are greatest, by 
employing logarithmic stretching: 

f 1 = log(i + e~^"' ) + Lx , (A6) 
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66.71 
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-3 


0.01 


43.48 



Table A4. The time taken for numerical dissipation to halt the 
simulation for a = 100 fm = 110, and a 50 x 100 grid, as a function 
of ggenl: xlrat, the geometric ratio in the rescaled coordinate r. 



= -log[l-(l-e-^'')y]. (A7) 

To implement these coordinates in ZEUS-3D, set 
ggenl :xlrat to {Xe^'' + l)(<3x-i)-'^ where < x < X and 
Lx controls the 'zoom'. Radial logarithmic scaling gives less 
grid resolution near the outer boundary whore the density 
is least and thus Atz is greater and ZEUS-3D runs for a 
longer time. 

When importing a Grad-Shafranov equilibrium 
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Figure A5. Isothermal atmosphere with dipolar magnetic field. 
The kinks in the magnetic field lines (solid) are due to numerical 
dissipation. Parameters: a = 50, fm = 60, f = 66, 150 X 150 grid 
{left); a = 100, fm = 110, f = 66,50 X 100 grid xlrat = 1.1 
(right). 

(PM04), whose grid is linear in y = coa9, into ZEUS-3D, 
whose grid is (in many cases) logarithmic in y, there is no 
trivial multiplicative factor relating yj and the required 
logarithmic 9 scaling through ggenl:x2rat. The problem 
can be overcome by rewriting the Grad-Shafranov code 
such that its grid is logarithmic in 6. This problem is an 
obstacle in certain sorts of numerical experiments e.g. the 
bootstrapping method in §4.2. 



